Gaussian resolutions for equilibrium density matrices 

Pavel Frantsuzov 1,2 , Arnold Neumaier 3 and Vladimir A. Mandelshtam 1 

1 Chemistry Department, University of California at Irvine, Irvine, CA 92697, USA 

2 Institute of Chemical Kinetics and Combustion RAS, Novosibirsk 630090, Russia 
3 Institut fur Mathematik, Universitat Wien Strudlhofgasse 4-, A-1090 Wien, Austria; 

WWW: http:/ /www. mat. univie. ac. at/^neum/ 
(February 1, 2008) 



m 

O 
O 

(N 
C 

00 



> : 

(N ■ 

vo : 
o . 

m ■ 
O ■ 

Oh! 

+-> ' 

G ■ 
cd 

3 : 



A Gaussian resolution method for the computation of equi- 
librium density matrices pr for a general multidimensional 
quantum problem is presented. The variational principle ap- 
plied to the "imaginary time" Schrodinger equation provides 
the equations of motion for Gaussians in a resolution of pr 
described by their width matrix, center and scale factor, all 
treated as dynamical variables. The method is computation- 
ally very inexpensive, has favorable scaling with the system 
size and is surprisingly accurate in a wide temperature range, 
even for cases involving quantum tunneling. Incorporation of 
symmetry constraints, such as reflection or particle statistics, 
is also discussed. 



Propagating Gaussian resolutions. Computing 
the equilibrium quantum density matrix 



p T := Z^e-P", 



Z := Tr e 



-0H 



(1) 



in dependence of the inverse temperature f3 — 1/kT, or 
finding an equilibrium property of the type 



(A) T ■= Tr p T A 



(2) 



for some operator A, is of great interest in statistical 
physics and molecular dynamics. The most accurate 
methods, based on a spectral resolution 



Pt = 



n)(n\ 



out to give a simple, very inexpensive and surprisingly 
accurate alternative to the other approaches. An exten- 
sion to noncquilibrium problems is possible and will be 
discussed elsewhere. 

Consider a ./V-body (or D — 3-ZV-dimensional) quantum 
system with Hamiltonian 



H = ~p r M- l p+ U(x). 



(3) 



Here p (x) define £>-dimensional column vectors of mo- 
mentum (coordinate) operators and M is the mass ma- 
trix. The potential function U (x) is assumed to be rep- 
resentable as a sum of terms, each being a product of 
a complex Gaussian times a polynomial. (With this as- 
sumption a commonly used plane wave representation of 
U(x) is covered.) 

Given a grid q n spanning the physically relevant region 
in configuration space we put 



(4) 



with (x\q n ) = 5{x — q n ). The resolution of identity can 
be approximated by 



1 = J d D q\q)(q\ ^^2w n \q n )(q n \, 



(5) 



where the sum is over all the grid points and the w n 
are quadrature weights defined by the inverse local den- 
sity of q n . (For large number of dimensions D a Monte 
Carlo procedure to generate q n may be adopted.) Writing 

e~P R — e~P H / Ie~P H / 2 and inserting Eq. 5, we obtain 



are extremely expensive and limited to very few degrees 
of freedom. Commonly used practical, albeit still very 
expensive, strategies involve path integrals [1]. An alter- 
native semiclassical method was recently proposed [2], 
but it is also expensive, while not very accurate. In this 
letter we modify the Gaussian propagation techniques of 
solving the time-dependent Scrodinger equation [3-6] for 
the present case, formally corresponding to propagation 
in imaginary time /3, and complement it by an impor- 
tant new idea: Monte Carlo-type sampling is replaced 
by the propagation of a Gaussian resolution of px- This 
brings our semiclassical technique formally closer to the 
full quantum treatment by a spectral resolution. It turns 



Wn(g„,/3/2|<2ri,/3/2), 



p T w Z 1 ^2w n \q n ^/ 2 ){q„,i3/2\ 



(6) 



(7) 



To find \q n ,j3/2) we note that it is the solution of the 
initial-value problem 



-H\q n .r), \q n ,o) = \g.n) (8) 
We solve Eq.(8) approximately using the ansatz [3-6] 

k,r) « |A(r)) (9) 
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with 

(x\X(t)) = cxp|7(r) 



-[x-q(T)] T G{T)[x-q(T)] 



(10) 



This is a Gaussian with center q, a real vector, width 
matrix G, a real symmetric and positive definite, and 
scale 7, a real constant, and A := (G, q, 7), a short-cut 
notation containing all the Gaussian parameters. 

With this ansatz the maximum number of parameters 
corresponding to the use of full DxD dimensional matrix 
G is (D + 2){D + l)/2. This number may possibly be 
reduced assuming weak coupling between certain degrees 
of freedom and setting the corresponding matrix elements 
of G to zero. The minimum number of 2D + 1 parameters 
would correspond to using a diagonal width matrix G. 

To solve for A = A(r) we follow the corresponding 
derivations [3-6] for the real-time dynamics of a Gaus- 
sian wavepacket by utilizing the variational principle 



dL 



= 



with the Lagrangian 



L=(x'(r) 



dr 



+ H 



Mr) 



Defining the two matrices 

K=(X'\X), H = (X'\H\X), 



Eq. 11 can be rewritten as 



d 2 K ■ dH 
^7 A + 



<9A<9A' 



dX' 



= 0, 



(11) 



(12) 



(13) 



(14) 



\'=X 



thus providing the equations of motion for the Gaussian 
parameters A = A(r). The use of the Gaussian wavepack- 
ets and the assumed form of the Hamiltonian (3) allows 
one to evaluate the corresponding matrix elements and 
their derivatives in Eq. 14 analytically. (Efficient nu- 
merical expressions will be published elsewhere, but see 
refs. [4,6]). Note that the commonly used linearization of 
the potential [3] reduces the complexity of the matrix el- 
ements evaluation, but also the accuracy, and is not used 
here. 

The delta function in the initial condition of Eq.(8) 
can be considered as a limit of the Gaussian (10) with 
infinite G. To avoid the singularities at r = we start 
instead at some small To where we approximate |g n ,T ) 
by a Gaussian that has finite width: 



(z|A(r )> 



' det M 

(27TT ) d 



(15) 



x exp 



^-{x-q n ) T M{x 



Qn) - U(q n )T 



To solve Eq. 14, we use the implicit integrator DASSL 
[7] , which has an error control and can be applied directly 
to Eq. 14. (Probably, a standard numerical integrator 
could be utilized, if Eq. 14 is explicitly solved for A.) 
Note that is small at large temperature kT = 1/0, and 
only a few integration steps are needed. As T decreases, 
the integration time becomes larger, and the variational 
approximation by Gaussians may become poorer. We did 
not encounter any special numerical difficulties, except in 
the cases when the Gaussian width was too small. 

Given (9), we can rewrite (7) as 

p T w Z- 1 Y,w n \X n (0/2))(X n (0/2)\, (16) 

and get for expectations 

(A) T ^ Z- 1 J2^n(X n (0/2)\A\X n (0/2)). (17) 

The evaluation of {A) requires the computation of many 
matrix elements of A between new Gaussian wavepackcts 
for every value of T . This can be done analytically if A 
is polynomial in p, x or the product of a polynomial and 
a complex Gaussian. 

Taking into account symmetry. Very often the sys- 
tem in question has symmetries, that is, the correspond- 
ing time-dependent Schrodinger equation conserves cer- 
tain symmetries. To show how to utilize this, we consider 
the example of reflection symmetry, Stp(x) = ip(—x), and 
assume that 



|g£ T > = ^(|gn,T>±S|g niT >) 



(18) 



is a solution of Eq. 8 at any r. Rewriting the resolution 
of identity, 



^2w n (\q+){q+\ + \q n ){q n \), 



(19) 



we obtain 

Z ^^2w n ((q+ 0/2 \q+ 0/2 ) + (q~ 0/2^-^2)), (20) 

P~ Z_1 S«'"(l< j 8/2>K9n,/J/2l + I^^K^/sl), 
(A) w Z- 1 Y J ^n{(qt N2 \A\q+ N2 ) + (q- (j/2 \A\q- 0/2 )). 

To approximate the solutions |g^ r ) we can replace the 
Gaussian in Eq. 9 by the symmetrized Gaussian: 



where 



|A ± ) = i=(|A>±£|A», 



S|G(r),g(r),7(r))=|G(r),-g(r), 7 (r)). 



(21) 



Therefore, the matrix elements needed in Eq. 14 can be 
evaluated by taking the appropriate linear combinations 
of matrix elements between single Gaussians. 
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Bose and Fermi statistics. Another important case 
of symmetry corresponds to the Bose or Fermi statistics. 
For a system of N indistinguishable particles in configu- 
ration space defined by the coordinates 



we should symmetrize or antisymmetrize the wavepackets 
according to the particle statistics. Let a = (aq, a at) 
define a permutation of particle indices, and denote by 
P a the permutation matrix with 



Exact quantum 
Gaussian propagation 
Classical 
H-K 




FIG. 1. Mean square displacement computed by four differ- 
ent methods for ID single well potential U(x) = ^x 2 + 0.1a; 4 . 
The semiclassical result labeled by "H-K" is taken from ref.l 
The difference between the present and exact quantum result 
is not seen in the graph. 



1.9 



I 



Exact quantum 
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FIG. 2. The mean square displacement for the symmetric 



double- well potential U(x) 
different methods. 



4 — Ax + x computed by four 



P(x3j — (^-Qi , • • ■ , Xa N ) 

A symmetrized wavepacket |A±) can be defined by 

(x|A±>=^sign ± (a)(P a x|A). (22) 

a 

Here boson statistics has sign_(a) = 1; in the case of 
fcrmion statistics, sign + (a) = 1 for even permutations 
and sign + (a) = —1 for odd permutations. The equations 
of motion (14) have the same general form. Since Eq. 22 
implies 

|A ± )=£sign ± (a)P*|A) 

a 

and one easily verifies 

P*\G,q,l) = |P*GP a ,P^, 7 >, 

here again, matrix elements with symmetrized Gaussians 
|A±) are simple linear combinations of suitable matrix el- 
ements with unsymmetrized Gaussians. Similar consid- 
erations apply to multiparticle systems with few indis- 
tinguishable particles, where the sum(22) has only a few 
terms and the calculations remain feasible. 

Numerical examples 
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FIG. 3. The density p(x) = (x\p\x) at kT = 0.5 computed 
by three different methods for the symmetric double-well po- 
tential U(x) as in Fig. 2. 

ID single-well problem. We first apply the method 
to a ID problem of ref. [2] with potential U(x) = 
i.x 2 + O.lrr 4 to compute the mean square displacement 
(x 2 ) — (x) 2 as a function of temperature T. A grid of 10 
equidistant points was taken in the interval — 5 < q n < 5. 
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(Here and in all the following examples the reported re- 
sults fully converged with respect to the grid size.) Each 
of the 10 Gaussians was then propagated by Eq. 14 start- 
ing with To = 0.01 up to r = 50. The result obtained by 
the present method is hardly distinguishable from the 
converged quantum calculation using diagonalization of 
the Hamiltonian in a large basis. In Fig. 1 these results 
are also compared to the classical Boltzmann average and 
to the semiclassical calculation using the Herman-Kluk 
propagator [2] . Note that in the latter case a much more 
expensive Monte Carlo method was used with 10 5 classi- 
cal trajectories, still resulting in relatively big statistical 
errors. 




kT 



FIG. 4. The mean square displacements computed by three 
different methods for the asymmetric 2D double- well potential 
U(x,y) = {x 2 - 2) 2 + 0.1a; + 0.5(y - 0.5a;) 2 + O.ly 4 

We note that a similar high accuracy was achieved for 
a 2D single well problem (not shown here). 

ID symmetric double-well problem. Here the 
method was applied to a problem with the potential 
U(x) = 4 — 4x 2 -I- x 4 . Now a grid of 14 points with 
— 3 < q n < 3 was used. The corresponding results using 



both the unsymmetrized and symmetrized gaussians are 
shown in Fig. 2. The agreement between the exact result 
(fully converged diagonalization of H in a large basis) and 
that computed by the present method is unexpectedly 
excellent even at quite low temperatures, well below the 
potential barrier. For the unsymmetrized case a signifi- 
cant deviation from the exact (and symmetrized) result 
occurs only at low temperatures where the small tun- 
neling splitting causes the quantum observables change 
rapidly with T. 

In Fig. 3 we show the density profile for the same 
system computed using the same set of 14 symmetrized 
Gaussians at kT = 0.5 together with the exact quantum 
and classical results. 

2D asymmetric double-well problem. To further 
demonstrate the method we apply it to a 2D problem 
with asymmetric double-well potential U(x,y) — (x 2 — 
2) 2 + 0.1a; + 0.5(y-0.5x) 2 + 0.1y 4 . An equidistant 16x16 
grid q n in a square box [—4; 4] x [—4; 4] was used. The 
results for the mean square displacements shown in Fig. 
4 are again surprisingly accurate, except for very low 
temperatures. 
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